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Abstract 

We have developed a coupled equations continuum model that explains the complex 
surface shapes observed in epitaxial regrowth on micron scale gratings. This model 
describes the dependence of the surface morphology on film thickness and growth 
temperature in terms of a few simple atomic scale processes including adatom diffu- 
sion, step-edge attachment and detachment, and a net downhill migration of surface 
adatoms. The continuum model reduces to the linear part of the Kardar-Parisi- 
Zhang equation with a flux dependent smoothing coefficient in the long wavelength 
limit. 
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1 Introduction 



The problem of the time evolution of the shape of crystal surfaces has a long 
history dating back to Mullins and Herring who considered relaxation during 
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annealing above the roughening temperature [1]. More recently, shape relax- 
ation below the roughening temperature has been studied extensively [3,4,5]. 
Below the roughening temperature the problem is complicated by the need to 
keep track of the dynamics of atomic steps and the fact that the surface free 
energy of crystal facets is singular. Biasiol et al. [6] have extended the theory 
of shape relaxation below the roughening temperature to include the effects 
of atom deposition, and use this theory to explain the self limiting V-grooves 
observed in organo-metallic chemical vapor deposition (OMCVD) growth on 
corrugated GaAs substrates. In this paper we present a new continuum model 
which we use to interpret measurements of the shape of corrugated GaAs 
(100) surfaces under growth conditions which do not produce faceting. Facets 
are not present in our experiments due to atomic scale roughness associated 
with atom deposition in the island growth mode, and the fact that the sur- 
face topography is sufficiently weak that the surface slope does not reach the 
low energy [111] facets. We show that this model reproduces the surface mor- 
phology that develops during molecular beam epitaxy (MBE) regrowth on ID 
surface gratings. 



2 Conventional Modeling of Weak Surface Texture 

The evolution of long wavelength surface structures during GaAs homoepitaxy 
can be described by the Kardar-Parisi-Zhang (KPZ) equation [1,2]: dh/dt — 
vV'^h + |(V/i)^ + F + ri{x, t). The coefficients in this equation are constants 
characterizing the microscopic atomic processes. The source term r]{x, t) simu- 
lates the random arrival of atoms at an average rate F. According to this equa- 
tion, a textured starting surface will develop parabolic mounds that smooth 
with time separated by V-shaped valleys. Recent experimental work [7] has 
shown that the KPZ equation provides an accurate description of the mor- 
phology of epitaxially grown GaAs layers for surfaces with local slopes < 3°. 
The agreement with this simple continuum model suggests that the anoma- 
lous effects associated with the singular free energy of crystal facets are not 
important for the growth conditions in question. 

In the case of GaAs molecular beam epitaxy (MBE) growth in which there 
is no re-evaporation, the simplest explanation for the linear term in the KPZ 
equation is that it is due to an inverse Ehrlich-Schwoebel (ES) effect [1] in 
which surface adatoms approaching a descending step are more likely to de- 
scend over the step rather than being reflected from it, due to a step edge 
potential barrier. This creates a downhill flux of adatoms (j oc —Vh) and 
a smoothing term {dh/dt oc —V • j) identical to the first term in the KPZ 
equation [1]. In practice, the atomic scale dynamics is complex with surface 
reconstructions, complicated step edge geometries, and a two component (Ga, 
As) surface [8,9,10]. The sign of the ES effect in GaAs is controversial [7,11,12], 
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Fig. 1. Light scattering during growth corresponding to surface power spectral den- 
sity at 41 /xm~^ showing the effect of atom deposition on the smoothing rate. The 
sample roughens during a temperature ramp to remove the surface oxide at about 
5 minutes in the figure, which is followed by relatively fast smoothing during a high 
temperature (620° C) anneal for about 7 minutes, and then slower smoothing during 
annealing at growth temperature (550°C). 

but we show below that a negative ES effect, favoring downhill flow (stable 
growth) is consistent with the experimental data. 

The nonlinear term in KPZ is associated with growth along the outward nor- 
mal, as in chemical vapor deposition. In this case, A should be equal to the 
growth rate F. However, the value for A needed to simulate the experimental 
results is almost two orders of magnitude larger than F [7]. Also, the KPZ 
nonlinearity is non-conservative, whereas MBE growth is conservative with a 
growth rate that is independent of the surface shape. 

In addition, the KPZ description with constant coefficients is not consistent 
with experiments which show that the smoothing rate depends on the growth 
rate. For example, in Fig. 1 we show the scattered light intensity from a GaAs 
surface during an interruption in growth on a randomly textured substrate. 
The intensity of scattered light is proportional to the power spectral density 
of the surface topography at a spatial frequency q determined by geometrical 
factors [7] (in this case q—Al //m~^, corresponding to a lateral surface length- 
scale of about 150 nm). For low amphtude surface textures, in the KPZ model 
the surface should smooth exponentially with a characteristic rate given by 
vq^ where q is the spatial frequency of the surface roughness [1]. As shown 
in the inset of Fig. 1, the smoothing rate responds immediately to changes 
in the growth flux; it is faster during deposition and slower during annealing. 
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Fig. 2. AFM images of (a) a sample quenched (fast cooled) after 69 minutes of 
growth at 600° C and (b) a sample annealed for 15 minutes at growth temperature 
595° C after 40 minutes of growth. 



suggesting that v is flux dependent. This continued smoothing of the surface 
in the absence of an atom flux indicates that the physical mechanisms at play 
on the surface still favor a net downhill migration of surface adatoms, even 
after the flux of atoms from the vapor has been turned off. 

Insight into why the smoothing rate depends on the flux can be obtained by 

comparing an atomic force microscope (AFM) image from a sample which is 
fast cooled (quenched) after growth with the AFM image of a surface which has 
been annealed (see Figs. 2a and b). The quenched sample (a) is covered with 
small islands, whereas the annealed sample (b) has broad terraces with few 
islands. The small islands must coalesce into the step edges during annealing. 
The kinetic barrier to the adatom coalescence into the step edges, causes the 
growth process to be non-local in space and time, in contrast to KPZ. A high 
density of steps at one location that absorb adatoms will affect the adatom 
density and hence the growth rate at another nearby location. 



3 Coupled Growth Equation Model 



The growth phenomena discussed above can be explained in a natural way if 
we extend the growth model to include the adatom dynamics explicitly with 
two coupled growth equations (CGE) [13]: 
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2Dn^ + {pDn - k)S. 



(lb) 



(la) 



Eqn. la is a continuity equation for the adatom density n with source and 
sink terms, while Eqn. lb describes the time dependence of the surface height 
h, which depends on the dimer nucleation rate and the net adatom attach- 
ment rate at steps. The constants are defined in atomic units as follows: F - 
deposition rate from the vapor, D - adatom diffusion coefficient, S - density 
of steps, K, - rate of thermal evaporation of atoms from step edges into the 
adatom phase, and < /3 < 1 is the sticking coefficient for an adatom at a 
step edge. An adatom is defined as a diffusing unit on the surface, which could 
be a Ga atom or a Ga-As complex. We also define: 



where in Eqn. 2a, j is the surface current of adatoms and < ^ < 1 is a 
proportionality constant which describes the net downhill drift of adatoms. 
The second term in Eqn. 2a represents adatom diffusion. In Eqn. lb, any 
adatom that attaches to a step edge is assumed to have incorporated into the 
film. The downhill drift parameter ( can be positive or negative: a positive 
value favors downhill drift of adatoms, consistent with the surface smoothing 
that is observed experimentally for GaAs (001) [11,7] (and also consistent with 
a negative value for the ES energy barrier) . 

In Equation 2b, we present a physically plausible hypothesis for the depen- 
dence of the rms step density on the surface slope. In this expression the 
random local surface slope associated with the growth-induced step density 
5*0 is added to the deterministic macroscopic surface slope Vh. Since the lo- 
cal slope associated with the background step density 5*0 is random, the two 
terms add in quadrature. In Equation 2b, we assume that the background 
step density is independent of the macroscopic surface slope. We expect 5*0 
to depend on temperature and deposition rate, and on time in the case of 
growth interrupts (see Fig. 2a and b) [14]. The simple picture of a surface 
consisting of fiat terraces separated by atomic steps, can be expected to pro- 
vide a good description as long as the surface slope docs not reach the next 
low index crystal planes, namely (110) and (111). These planes are 45° and 
54.7° from the surface normal, and beyond the range of surface slopes that we 
have explored experimentally (< 30°). We assume that the density of random 
steps 5*0 is independent of the topography. In this approximation, the average 
step density is proportional to the rms value of the local surface slope. The 
expression for S is then constructed by averaging over the random orientation 
of the local slope, and the rms step density is given by the incoherent sum of 



j^-D{CnVh + Vn), 



(2a) 
(2b) 
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Fig. 3. Film thickness dependence: (A) AFM scan lines for regrowth on 100 nm 
deep gratings oriented perpendicular to the [110] direction; (B) Scan lines from 
CGE calculation; (C) Scan lines from 2D kMC simulation of 10 nm high grating 
structure, where one At equals 5.6 ML of growth. All offsets arbitrary. 



the two contributing factors. 

For low amplitudes and long wavelength (V/i < Sq), the adatom density will 
be nearly constant as a function of position and time, and approximately equal 
to riQ = {F + uSq)/ I3DSq. In this case, Eqns. la, lb reduce to. 



dh 
'dt 



(3 \So 



(3) 



This reproduces the linear part of the KPZ equation and shows explicitly the 
dependence of the linear smoothing coefficient v on the deposition rate and the 
downhill drift parameter C,. In addition, it shows that in the absence of growth 
{F = 0) the linear smoothing term is independent of the background step 
density 5*0. This agrees with the light scattering data in Fig. 1, which shows 
that the smoothing rate is relatively constant during a growth interruption 
even though the AFM images in Fig. 2 indicate that the step density drops 
dramatically during annealing. Extending Eqn. 3 to higher order, one finds 
non-linear terms with higher order spatial derivatives. We speculate that the 
higher order nonlinear terms can be approximated by the KPZ nonlinearity 
over a limited spatial frequency range if the surface topography is not too 
large. 
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4 Textured Surfaces: Film Thickness Evolution 

Growth on substrates with larger amplitude surface slopes, up to ~ 30°, show 
complex surface shapes before evolving into parabolic mounds, as shown in 
Fig. 3a. At intermediate times the valleys are V-shaped with concave rather 
than convex sidewalls and distinct shoulders near the top of the sidewalls. 
Note the absence of (100) facets which are predicted theoretically for annealing 
below the roughening temperature in the absence of deposition [5] . Equations 1 
and 2 can be solved in seconds with a finite difference scheme and a coupled 
differential-algebraic system solver, and a ID solution is shown in Fig. 3b with 
parameters adjusted to match the experimental data in Fig. 3a (see Table 1 for 
parameters). The agreement with the experimental surface shapes is striking. 
In particular, the model reproduces the inverted "Gothic window" shape of the 
valley for the 600 nm growth and the KPZ-like cusps in the 2600 nm growth 
where the grating amplitude has reduced sufficiently so that the structure is 
described by the KPZ equation. 

A continuum model cannot include the microscopic details of the atomic scale 
phenomena, such as the geometry and density of step edges. We therefore com- 
pare the continuum model in Eqns. 1 and 2 with a kinetic Monte Carlo (kMC) 
simulation, which includes the same physical processes that are included in 
the Eqns. 1. We use a 2D, cubic grid, one-component, restricted solid-on-solid 
(SOS) model, with nearest-neighbor interaction. Each atom bonds to the sur- 
face with an activation energy Eact = Egub + n^Eiat, where m is the number 
of lateral neighbors [15]. The kMC simulations produce a random step distri- 
bution automatically due to the statistical nature of the model. In kMC, the 
binding energy for an atom at a step edge depends on how many neighbors it 
has (~ mEiat), whereas in the CGE continuum model a single average value 
is used for the step edge binding energy. 
Table 1 

Parameter table for CGE calculations. Atomic units were used with a lattice con- 
stant of 0.3 nm was used. 
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SOS simulations of MBE growth by kMC are limited by available comput- 
ing power to small scale structures, and become intractable for realistic, high 
temperature growth scenarios where 2D systems have sides up to microns and 
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growth times on the order of hours. In Fig. 3c, we show a kMC simulation for 
a surface grating that is somewhat smaher than the experimental structures. 
The simulated grating profiles in Fig. 3c were obtained by projecting 2D kMC 
simulations onto a line at each time point by taking the average elevation 
perpendicular to the scan hue. In this simulation, Egub—^-'^^ eV, Eiat— 0.35 eV 
and an ES step-edge barrier of Ees=-0.05 eV was used for the downhill drift 
mechanism. The agreement with the experimental shapes is excellent, repro- 
ducing all of the main features, except they are on a smaller size scale. The 
substrate and lateral binding energies are similar to values reported earlier 
in the interpretation of RHEED data [16,17] and compatible with the fitting 
parameters found in the continuum model. It is plaiisible that similar shapes 
could be obtained for the larger size scales relevant to the experiments by scal- 
ing the parameters appropriately. In the case of the CGE model (Eqns. 1, 2) we 
find that the parameters can indeed be scaled to yield similar surface shapes 
at different length scales [18]. 



5 Textured Surfaces: Temperature Evolution 

In Fig. 4a, we show the dependence of the surface topography on growth 
temperature, for a fixed layer thickness together with (b) the simulated surface 
topography using Eqns. 1 and 2 and parameters as outlined in Table 1. The 
experimental data is obtained from growths on 100 nm deep gratings oriented 
perpendicular to [110]. This is the faster diffusion direction in this material 
system [7], and depends on the As2/Ga ratio during growth, which was equal 
to three for all growths discussed in this work. This observation is consistent 
with the values used for the downhill drift parameter in our calculations, 
where the best fits were obtained using a larger ^ when the gratings were 
oriented perpendicular to [110] (Fig. 4b) than for gratings perpendicular to 
the [110] direction (Fig. 3b). The diffusion constant D, however, was considered 
isotropic in all calculations in this paper. There is some uncontrolled variation 
in the pitch and depth of the gratings in the experimental data in Fig. 4a. 
Nevertheless, the CGE model reproduces the main features in the temperature 
dependent data, namely the small secondary mound in the valley at 500°C, the 
KPZ-hke cusp at 550°C and the inverse Gothic window shape for the valleys at 
610°C. The shoulders at the edges of the ridges at 610°C are also reproduced 
by the model, although they are not as sharp as in the experimental data. 

The parameters used in these calculations are based on the same energies used 
in the kMC simulations in Fig. 3c. The diffusion constant is related to the sub- 
strate binding energy through the expression D = {2kT/h) exp {—Esub/kT). 
The step edge detachment rate is calculated from k, = D exp{—rhEiat/kT), 
where rh is an average number of neighbors for atoms at a step edge which we 
set equal to 2.25. The dechning value of Sq with temperature is reasonable; 
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Fig. 4. Growth temperature dependence: (A) AFM scan lines for regrowth on 100 
nm deep gratings oriented perpendicular to the [110] direction; (B) CGE calculation. 
The grating pitch is 5/xm. All growths are 1 hour at 0.8ML/s. All offsets arbitrary. 



one might also expect ( to decrease weakly with temperature. Satisfactory fits 
to the data were also obtained with a larger activation energy for D (1.8 eV 
rather than 1.25 eV) together with a smaller prefactor and somewhat different 
(yet still physically reasonable) temperature dependences for the other param- 
eters. Experimental and theoretical work suggests that the activation energy 
for D is in the 1.5-2.0 eV range. [16,17,19,20]. 



6 Conclusions 



We have shown that the complex surface morphology which develops during 
epitaxial regrowth on patterned GaAs (100) substrates, can be explained by 
the dynamics of the deposited adatoms, including step edge attachment and 
detachment, adatom diffusion, and downhill drift. Although we attribute the 
downhill drift to a negative Ehrlich-Schwoebel barrier we cannot rule out the 
possibility that this effect is caused by some other mechanism. This analysis is 
specific to GaAs, but the concepts are generic and may be applicable to other 
systems. 
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